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We study the evolution of density in an expanding Bose-Einstein condensate ttiat initially has a spatially 
varying phase, concentrating on behaviour when these phase variations are large. In this regime large density 
fluctuations develop during expansion. Maxima have a characteristic density that diverges with the amplitude 
of phase variations and their formation is analogous to that of caustics in geometrical optics. We analyse in 
detail caustic formation in a quasi-one dimensional condensate, which before expansion is subject to a periodic 
or random optical potential, and we discuss the equivalent problem for a quasi-two dimensional system. We 
also examine the influence of many-body correlations in the initial state on caustic formation for a Bose gas 
expanding from a strictly one-dimensional trap. In additon, we study a similar arrangement for non-interacting 
fermions, showing that Fermi surface discontinuities in the momentum distribution give rise in that case to sharp 
peaks in the spatial derivative of the density. We discuss recent experiments and argue that fringes reported in 
time of flight images by Chen and co-workers [Phys. Rev. A 77, 033632 (2008)] are an example of caustic 
formation. 
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I. INTRODUCTION 

Focussing of rays and the associated phenomenon of caus- 
tic formation are both well known in optics 1 1]. For the phase 
screen model caustics have been studied extensively by Berry 

• In this model a monochromatic plane wave encounters a 
thin screen, located on the plane z—Q and having coordinates 
C and rj within the plane. The screen impresses on the wave 
a phase d{C, rj), which may be deterministic or random. In 
either case, for strong variation of this phase, a wave propa- 
gating in the z-direction and passing through the screen will 
develop large intensity variations. Specifically, observation 
of the wave intensity at a point sufficiently far beyond the 
screen will reveal a pattern of bright lines. These are caus- 
tics. Within geometrical optics, light intensity on caustics di- 
verges. Diffraction effects smooth these singularities and dec- 
orate caustics with interference fringes. 

A similar phenomenon can occur with matter waves associ- 
ated with propagating clouds of cold atoms, especially if these 
form a Bose-Einstein condensate (BEC). The purpose of this 
paper is to develop a theory of caustics for cold atoms and to 
discuss the experimental conditions for their observation. The 
arrangement we consider differs substantially from that in op- 
tics. Caustics develop during the expansion of an atomic cloud 
released from a trap. The corresponding matter wave is not at 
all monochromatic, and time assumes the role of the spatial 
axis of propagation in optics. The mechanism of impression 
of the phase is also different. One possibility is to create den- 
sity variation in the trap. During the expansion this initial 
density modulation, in combination with strong non-linearity, 
produces a space dependent velocity field, thus impressing a 
phase on the BEC In Sec. |II]we shall discuss further this 
mechanism, and the way it can lead to formation of caustics, 
for quasi-one dimensional systems. Another possibility for 
impressing phase variation is by applying to a trapped conden- 
sate a short pulse of a space-dependent potential. Immediately 



after the pulse, the trapping potential is switched off so that 
the condensate starts its free expansion with the phase vari- 
ation generated by the pulse. We employ this mechanism in 
Sec. Uni where we discuss the strictly one-dimensional case. 

In both examples, the geometrical optics limit is the regime 
in which the impressed phase has spatial variations that are 
much larger than unity. Characterising these phase variations 
by an amplitude 6*0 and a spatial scale Rq, our central con- 
clusion is that for 6*0 > 1 an expanding condensate devel- 
ops caustics after an expansion time of order t* = mR^/Wo, 
where m is the atomic mass. The density close to caustics di- 
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verges with 6*0, as 9^ in the simplest case, and may therefore 
be much larger than it is in the background between caustics. 

The development of density fluctuations as a consequence 
of initial phase fluctuations has been studied in previous ex- 
perimental and theoretical work, both for the case when these 
phase fluctuations are thermal in origin |3], and for the case 
in which they arise from a disordered potential applied to the 
trap|4]. In this earlier work, however, the significance of 
was not identified, and the method used were applicable for 
6*0 <C 1, when density fluctuations never become large. 

In outline, the organisation and main results of this paper 
are as follows. In Sec. we treat quasi-one dimensional 
systems with an initial density modulation, using the Gross- 
Pitaevskii equation to describe the conversion of density to 
phase modulation. Relative density variations generated in a 
trap by a potential with amplitude Vo small if Vb is much 
smaller than /i, the chemical potential. Nevertheless, they 
may lead to phase fluctuations that are large, since with ra- 
dial frequency ll!±_ we find 6*0 ^ Vo/^^-L- We suggest that 
the large density contrast measured recently jsl] for a conden- 
sate expanding from a trap with a disordered optical poten- 
tial should be understood as an example of caustic formation. 
We also examine behaviour as Vq is increased to values larger 
than /i, inducing fragmentation of the condensate. We show 
that caustics are not formed in the expansion of a highly frag- 
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mented condensate, but find that the value of Vq above which 
they are eliminated is larger than the threshold for fragmenta- 
tion. In Sec. inilwe examine the effect of many -body corre- 
lations in the initial wavefunction on caustic formation, tak- 
ing these correlations from the Lieb-Liniger model. We find 
that behaviour is controlled by the value of the healing length 
^: taking Oq ^ 1, caustics survive if interactions are weak 
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and 9q' $/i?o ^ 1, while in the opposite limit they are sup- 
pressed. We also consider, in Sec. IIVI expansion of a conden- 
sate from a two-dimensional trap in which there is a smooth, 
spatially varying potential. In this geometry caustics form a 
network of intersecting curves for huj± ^ Vq <C /i. Segments 
of these curves are eliminated for Vq ^ /i, although with- 
out any sharp signature of the percolation transition that takes 
place at fi — fic for the condensate in the trap. For ^ /ic 
caustic formation is suppressed altogether, as in quasi-one di- 
mension. We discuss briefly behaviour for fermion systems in 
Sec.|V] and close with a summary in Sec. IVIl 



II. QUASI-ONE DIMENSION: MEAN FIELD THEORY 

In this section we consider a strongly anisotropic BEC, ini- 
tially confined in the radial direction by a harmonic trap with 
frequency uj±. In experiments there is also weak confinement 
in the axial direction, with frequency uj^ ^ uj±- The axial 
confinement will be neglected, which implies that our consid- 
erations are limited to times i <C 1 /uJz after the release of the 
condensate from the trap. We use z and p as axial and radial 
coordinates. 

Prior to its release the condensate is in its ground state (we 
assume zero temperature) in the presence of a radial confin- 
ing potential imcj^p^ and a z-dependent potential V{z). In 
order to clarify the mechanism of caustic formation we shall 
start by treating a periodic potential V{z) — Vo cos kgz, and 
then proceed to the case of a disordered potential. In the lat- 
ter example we use Vq to denote the characteristic amplitude 
of the random potential, and Rq its correlation length, with 
i?o ^ 1/fco for the two cases to be comparable. We assume 
a smoothly var ying poten tial, in the sense that koa± <C 1, 
where ax = y^2fi/mLu''^ is the radius of the BEC in the trap, 
with fi the chemical potential, assumed much larger than liuj±. 

The potential V{z) produces density modulations of the 
BEC in the trap. Within the Thomas-Fermi approximation 
the ground state density is 



naip, z) = - [fi- V{z) - \muj\p^ 
9 \ 2 



(1) 



where g is the coupling constant for the non-linear term in the 
Gross-Pitaevskii equation. 

At time t—0 all potentials are switched off and the conden- 
sate expands according to the equations L6J 
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div nv = 
1 



at 

TO— + V ( —mv^ + gn ] =0 



(2) 
(3) 



where v = {vp, Vz) is the condensate velocity, related to its 
phase 9hy V = :^V6'. These equations are to be solved with 
the initial condition Eq. ([T]i, supplemented by the requirement 
that the velocity field v = 0, or equivalently that the phase is 
uniform, at the start of the expansion. 

The condensate undergoes rapid radial expansion, accord- 
ing to the standard scaling picture |6], but due to the initial 
density modulation it also develops an axial velocity compo- 
nent Vz{z, t). The latter is governed by the z-component of 
Eq. (O . During the initial stage of radial expansion we may 
neglect the kinetic energy ^mv"^ compared to the interaction 
energy gn, and obtain ||4|] 



Vziz,t) 



1 dV{z) 



This corresponds to an impressed phase 
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-V{z) 
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at any time lying in the window ^ io ^ l/w^. Both 
this phase and the consequences of non-zero Vz{z, tp) for the 
subsequent time evolution have been discussed in Ref. 
but the theory developed there applies only to systems with 
sufficiently weak initial density modulations. In the follow- 
ing we show that the typical phase magnitude 9q ^ Vq/ huj± 
is the relevant parameter: for 9o 1 only weak density 
modulations develop at later times, but for ^ 1 density 
modulations arise that are comparable to the average density, 
while in the limit 6*0 S> 1 caustics appear Large effects 
of this kind have been observed in a recent experiment iH. 
As we shall see, they develop at characteristic times of order 
t* = m/hkQOo which is parametrically larger than l/uj± by 
the factor /i/Vofega'^. 

During the second stage of expansion, at times much larger 
than to, the nonlinearity of the Gross-Pitaevskii equation can 
be neglected and we arrive at the problem of linear time evo- 
lution of the BEC wavefunction with, as an initial condition, 
the impressed phase 6{z) given by Eq. (|5]l. The wavefunction 
can be factorised into radial and axial parts, as 



^ip,z,t)^^pA)^{z,t) 



and the density is 



nip,z) = mp,t)\'\^iz,t)\' 



(6) 



(7) 



The behaviour of |<i>(p, <)p is well established and given 
by the radial scaling function of Ref. j^. Our concern in 
the following is with the axial part, \tp{z,t)\'^, which gives 
the density at point z and time t, normalised by the radial 
factor |$(p, The function ip{z,t) satisfies the linear 

Schrodinger equation 



ih 



~dt 



with the initial condition 



ip{z,to) 



2m dz^ 



(8) 



(9) 



3 



This form for ip{z,to) neglects the initial modulation of the 
BEC density, which is justified if Vq ^ /x so that the relative 
variation at time t — to is much smaller than unity. Since 
we are interested in density variations, at later times, that are 
of order unity or larger, we neglect these small initial modula- 
tions, except for their crucial part in producing the phase 9{z). 
For most of this section we impose the condition Vq <^ fi (but 
not the more stringent one Vq ^ huj± implicitly assumed in 
Ref. [4]). The opposite case Vq > fi corresponds to a frag- 
mented BEC and will be treated separately towards the end of 
the section. 

The solution to Eqns ([8]l and (|9]l for i ^ to is 



dz' exp 



2ht 



{z - z'f + i0{z') 



(10) 

where the impressed phase 9{z) can be an arbitrary function 
of z. 

Let us now consider 9{z) — 6'ocosfco2 and introduce the 
dimensionless variables z — k^z and i = t/t* = hkoOot/ni, 
so that 



dC e^p[i9o^iC,i,i)] , (11) 



with 



(p{C,z,i) 



2t 



cosC 



(12) 



The relative density \'(l){z,t)\'^, initially unity, acquires spa- 
tial variations with the passage of time. For <C 1 rela- 
tive density variations remain small for all times, as is clear 
from expanding the factor cxp(i6'o cos C) in the integrand of 
Eq ( fTTT ). This is the regime considered for a random potential 
in Ref. |4]. We study the opposite case, 6*0 ^ 1, when the rel- 
ative density can develop large modulations and caustics can 
be formed. This is the regime encountered in the experiment 
of Ref. |5]. 

We will be interested in the form of the relative density at 
a given instant t. This gives the z-dependence of the density 
of an expanding atomic cloud, measured with a probe beam 
perpendicular to the condensate axis. For t ^ 1, relative den- 
sity modulations are small. They grow linearly with t and are 
oscillatory in z with period 2ti. Growth of density maxima 
with time culminates, for 6q ^ 1, in the formation of caustics 
at times t > 1. In this regime Eq. ( fTTT i can be evaluated by 
the stationary phase method. Caustics are determined by the 
vanishing of not only the first but also the second derivative of 
the phase Lp{C,, z, i) with respect to C,: 



dC i 



- sin C = 



cos C = 



(13) 



(14) 



These equations have a simplegeometric interpretation, in 
complete analogy with optics lUI?!], which is illustrated in 
Fig. [T] Eq ( fTST l defines rays of atoms, in the sense that 






FIG. 1: Paths z{t) followed within the geometrical optics approxi- 
mation by atoms expanding from a condensate with an initial phase 

6[z) = 6o cos{z). 



atoms emerging from a point C have a velocity — sin C, given 
by the derivative of the phase cos C,, and so reach the point 
z = C,—t sin at time t. The roots of Eq. ( fT4l ) identify singu- 
lar points C,n- Vanishing of d'^ip{C, z, t) at these points means 
that the emerging rays are focussed on the observation points. 
Thus, for a given i, one will observe caustics, as bright spots, 
at points z„ — C„ — isin^n- Fort 1, z„ « ^+n7r— (— l)"t 
with n = , ±1 , . . .. The density at these points is found 
by computing the integral in Eq. (fTTT) . Within the stationary 
phase method, since at the singular points z — Zn the first two 
derivatives of z,i) vanish, the integral is controlled by the 
third derivative, d^tp{(^, z, t), and has a value proportional to 

9q resulting in a large relative density |-0(z,t)p ~ ^'o^'^- 
For times t 3> 1 the density in the vicinity of caustics decays 
as t^^, as is clear from the prefactor to the integral in Eq. ( fTOl ). 

At caustics one can observe spectacular diffraction effects 
[IS]. On slight deviation Sz ^ 9q from the point z„ the den- 
sity displays a sharp drop, followed by aperiodic oscillations. 
The detailed shape of these oscillations can be found by fur- 
ther studying the integral in Eq. ( fTTT ). We shall not do this here, 
but rather remark on the picture between caustics, at a point z 
well separated from all z„. In this case only Eq ( fTsT l remains 
to be satisfied, and since it can have several solutions, several 
saddle points will contribute to the integral in Eq. ( fTTT i. Each 

— 1/2 

saddle point contributes a term of order 9^ ' which cancels 
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the prefactor of 9q . Interference between different contribu- 
tions — that is, between different rays arriving at the point z 
at time i — results in a density pattern with variations of order 
unity and spacing between fringes of order Az ^ 1/9q. 

Our discussion has been limited to caustics of the simplest 
kind, referred to as folds in Ref. [2]. More singular caustics, 
known as cusps can also occur These stem from rays emerg- 
ing from the points (n = 2?™ with n = , ±1 , . . . and are 
visible only at a patricular time instant, i — 1. The point is 
that for (n — 2?™ and t = 1 not only the first two but also 
the third derivative, dQtp{C, z, t) — sinC„, vanishes. The in- 
tegral in Eq. ( fTTT i is then controlled by the fourth derivative. 



d*(p{(^, z, t) = 1, and is of order 



-1/4 



which implies that 
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FIG. 2: Relative density |?/>(2, t)p as a function of position i, for 
Oo = 30 and i = 0.5 (top), 1 (middle), and 2 (bottom). 
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FIG. 3: Relative density |?/>(i, t)p as a function of position z, for 
t = 3 and So = 10 (top), 30 (middle), and 100 (bottom). 



the reduced density at the points z„ = 27rn is at this time 
~ ^y^- Studying the integral in Eq. ( fTTT i in more 
detail, for values of t and z near t — \ and z — 2?™ one can 
identify lines in the z-t plane on which t)p drops from 
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its maximum value, of order 9^ , to values of order 9^ . 
These lines have cusps at the most singular points, located at 
z — 2?™ and t = 1. 

We illustrate these ideas using a numerical evaluation of 
Eq. (fTTI) to obtain the relative density The results 

are shown for several values of i at fixed in Fig. 2, and for 
several values of 6*0 at fixed t in Fig. 3. In the large 6*0 limit, 
caustics are located at: z = for i — 1 al\ z ^ ±0.685 
for < = 2; and ?A z — ±1.598 for t — 2>. The two figures 
show that \'4){z, i)p has peaks near these points, which grow 
in amplitude and become increasingly well defined with larger 
9o- 

We now turn to a disordered potential V{z), which we take 
to be Gaussian distributed with zero mean, amplitude Vb, cor- 



relation length i?o, and correlation function 



V{z)V{z') = Vo'f 



Rn 



(15) 



Here and elsewhere, we use an overbar ... to denote a dis- 
order average. The function f{z) has unit amplitude and unit 
range, so that /(O) = 1 and f(z)^0 for z ^ 1. This implies 
for the impressed phase that its mean is zero and 



e{z)9{z') = e'j 



Rn 



(16) 



with 9q ~ TTVo/2hLLjj_. The experimentally realised disor- 
dered potentials for cold atoms are optical speckle potentials 
(see Ref. ^ for an extensive discussion) whose correlation 
function is 



fiz) 



sin^(z) 



(17) 
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Eq. ( fTOb for iIj{z, t) applies also for a disordered potential and 
the mechanism for caustic formation is qualitatively the same 
as for the harmonic potential discussed above. Caustics are 
due to rays emerging from singular points z'^ on which there 
are zeros of two or more derivatives with respect to z' of the 
overall phase in Eq. (fTOl l. 



m 

mi 



iz-z'f + e{z') 



Since the impressed phase d{z') is now a random function, 
the formation of caustics of different types (fold, cusp or 
more singular) is a matter of probability. The typical time 
for caustic formation is t* — mR^/hOQ. Using the relations 
^mijj\a\ ~ fj, and 9o — TrVo/2hu!±, the formation time can 
be written in terms of the experimentally controlled parame- 
ters as 



Truj± Vb \ a 



It is quite straightforward to modify the theory of optical caus- 
tics [2, 7] to the case of caustics in an expanding BEC, but we 
do not pursue this further here. 

So far our discussion of caustics for the random case was 
qualitative and pertained to a specific, typical realisation of 
the disordered potential. We next present some quantitative 
analytic results for the second moment of the reduced density 
fluctuations. 



\^P{z,t)\^ = Sit) + l . 

The averaging restores translational invariance and so elim- 
inates the z-dependence. We denote the resulting (time- 
dependent) quantity by S{t) + 1 in order to emphasise the 
analogy with speckle patterns in optics. There S is called 
the scintillation index and it is a measure of spatial intensity 
fluctuations in the speckle pattern. For a uniform intensity S 
is clearly zero, and for the standard speckle pattern, created 
by great number of interfering waves with random phasesjn 
which the intensity has a Rayleigh distribution, S is unity f^. 
For the random screen problem, S was studied as a function 
of distance D from the screen in Ref . lIToll . It was shown there 
that, at first, S increases with D and reaches a maximum value 
larger that 1 but then, with further increase in D, the value of 
S drops and approaches 1 for large D. We briefly outline a 
similar calculation for our problem. 

Starting from Eq. ( fTOb we write |-!/)(z, as the product 
of four intergals. Performing the standard Gaussian average 
of the expression exp[i(6'(zi) — 0{z2) + dizz) — d(z4))] then 
yields 



S{t) + 1 = -^^ / / dxdye-"'°^y/'e-'o^'^^-y'> 



(18) 



The function K{x,y) is defined by 

K{i, v) = 2- 2/(5) - 2/(y) + f(i + v) + f{i - y) , (19) 
where f{x) is the correlation function introduced in Eq. ( fTSl ). 



For 00 <^ 1, exp{—9QK{x, y)) may be expanded as 1 
—9QK{x,y). Then S{t) is small, given by 



S{t) = 2e\ 



'271 



ds (sin + cos s'^)f{s 



(20) 

and saturates at the value 29q for t ^ Oq. This small Oq regime 
was studied in Ref. We turn to the more interesting case, 
Oq ^ 1, when caustics and related large interference effects 
can be observed. For this case no expansion in Oq is possible 
and the analytical treatment becomes tedious 1 10]. It is pos- 
sible, however, to identify three distinct contributions to the 
integral in Eq. ( fTSl l. each of which dominates at the appropri- 
ate time, and thus to derive an approximate expression: 



S{t) + 1 



— sin(6los/i)e-5T^o" 



'—7 / ds ssm{s/i)\n{s/9o)e 
Jo 



+2 



1 — erf 



/ 1 



\2^i 



(21) 



where /3 = |5|/(x)|i=o| and 7 = |9|/(i)|£=o|- For the 
speckle pattern these numbers are (3 — 1/3 and 7 = 2/45. 
The three terms in Eq. ( |2TI ) make the dominant contribution 
to the scintillation index at different times. For short times, 
t ^ 1, the first term dominates and approaches 1 as i — > 0, so 
that S{t) drops to zero. The last term dominates in the oppo- 
site limit, t ^ 1. It approaches 2, thus yielding the expected 
saturation value S{t) — 1. The most interesting term, how- 
ever, is the second one. It dominates for intermediate times, 
when t ^ 1, and it is proportional to In^^o. signalling the ap- 
pearance at this time of caustics and the associated large den- 
sity fluctuations. 

Let us emphasize that the two basic requirements for our 
treatment are huj± <C /i and Rq ^ a±. The first inequality 
ensures that the size and energy of the condensate, while in 
equilibrium in the trap, is dominated by the interactions i.e. 
by the nonlinear term in the Gross-Pitaevskii equation. The 
second inequality is required for the validity of the two-stage- 
expansion scenario as well as of the ray picture on which the 
physics of caustics rests. To this point our treatment of den- 
sity modulations in a BEC after expansion has been restricted 
to systems in which the relative amplitude of initial density 
fluctuations is small, which is the case when these are pro- 
duced by a potential with amplitude Vq <^ ^. It is in fact 
straightforward to generalise the discussion, to allow for an 
arbitrary value of Vq. A potential with amplitude Vq > n 
automatically implies impressed phase fluctuations with am- 
plitude do ^ 1. That in turn justifies a stationary phase treat- 
ment of caustic formation, and according to this treatment, 
the atoms that form caustics come from short segments of 



the condensate, which have width proportional to 



-1/3 



be- 



fore the second stage of the expansion. Moreover, the density 
in these segments remains approximately constant during the 
first stage of expansion, provided Rq ^ a±. To obtain the fi- 
nal relative density after expansion under these conditions, the 
value of \ijj{z, i)p calculated from Eq. ( fTOl i should simply be 
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FIG. 4: Schematic illustration of the effect of condensate fragmen- 
tation before expansion on caustic formation after expansion. Mid- 
dle panel: potential V{z) and chemical potential /i before expan- 
sion. Lower panel: condensate density li/jp before expansion. Upper 
panel: position z{t) reached by atoms after expansion as a function 
of starting point z, for three different expansion times, t = 0,ti and 
t2, with < ti < t2. Caustics appear at turning points of z{t) vs 
z, marked here for t2 at Za and zt- For the value of /i illustrated, the 
initial density is zero at points from which caustics would otherwise 
develop. 

multiplied by the initial relative density 0)p at a point 

z' determined from 

z' = z^ —d,,e{z') . (22) 

TO 

The consequences of this are most significant for Vq > /i, 
when the condensate before expansion is fragmented, having 
initial relative density \^p{z' , 0) ^ = for some values of z'. In 
particular, it turns out that caustics are completely suppressed 
for Vb ^ /i, because in this limit the initial density is zero at 
all points z' for which there is focussing of the emerging rays. 
To show this, note that focussing at time t of atoms from z' 
occurs if 

dl0{z') = -^^, (23) 

that is to say, d'^,9{z') must be negative. On the other hand, 
for Vq ^ /i, the condensate before expansion occupies the 
neighbourhood of minima of the potential V{x). In these re- 
gions d'l,V{z') and hence dl,0{z') are positive, while at the 
points where d'^,9{z') is negative, the initial density is zero. 
These ideas are illustrated schematically in Fig.lD 

We next comment on a recent experiment |5] which appears 
to satisfy the conditions for caustic formation. In this experi- 
ment ^/fiuj± — 5.6, Rq = 15/im, a± < 10/im and the largest 



density variations were observed for Vq = 0.5/i, in a time of 
flight image at ^toF — 8ms, which is significantly larger than 
l/uj±. This value for Vq corresponds to a phase amplitude 
^0 = 4.4 which gives the caustic formation time t* = 5ms. 
This is smaller than but comparable to the value of ItoF, and 
so the large density variations in Fig. 4(h) of Ref. ^ can most 
likely be attributed to caustics. While the value 9o — 4.4 does 
not lie deep within the large Oq regime, the experiment was 
not designed for caustic observation. The results presented 
above should facilitate the optimisation of such experiments. 
We note further that at Vq = fi, the largest potential ampli- 
tude for which results are reported in Ref. |5], the condensate 
before expansion appears to be fragmented (Fig. 4(i) of [5]) 
while there are no large density fluctuations after expansion 
(Fig. 4(i) of [5]). Both features are consistent with the sce- 
nario presented in our discussion of Fig. |4] 

A related experiment has been described in Ref. 101 to- 
gether with a theoretical discussion appropriate for small Oq. 
Some of the parameters for this experiment take similar val- 
ues to those of Ref. 11: iJ,/huj± = 6.8 and Vq = 0.41^, 
giving Oq = 4.3. However, the length scales Rq = 1.7/im and 
aj^ — 1.5^m are about an order of magnitude smaller than 
the corresponding ones in |5]. These lengths are also much 
smaller than the reported resolution (8.5/im) of the imaging 
system and of the density fluctuations after expansion that are 
illustrated in Fig. 1 of Ref. |14j]. We note in addition that the 
time of flight used in 0] (in the range 5- 17ms) is about an 
order of magnitude larger than the value t* ~ 1ms we cal- 
culate from experimental parameters. Therefore caustics in 
this experiment should form and decay before time of flight 
images are made, and should have a smaller spacing than the 
resolution of the imaging system. For these reasons we do not 
expect the theory we have presented to apply directly to the 
measurements of Ref. iQl. 



III. STRICTLY ONE DIMENSION: MANY BODY 
CORRELATIONS 

In this section we discuss the possibility of caustic forma- 
tion, in a system of interacting bosons, beyond the mean field 
approach. We assume here that huj± is much larger than the 
characteristic interaction energy so that, with respect to trans- 
verse motion, all N atoms in the trap reside in the ground 
state xq [p) of the harmonic oscillator, forming a strictly one- 
dimensional system. There is no confining potential in the ax- 
ial direction, in which particles move on an interval L with pe- 
riodic boundary conditions. The axial motion is controlled by 
an effective one-dimensional Hamiltonian - the Lieb-Liniger 
model ifm [T2I1 . The chemical potential of the system, in 
equilibrium in the trap, differs from hjj±^ only by a small, in- 
teraction induced correction. This correction is crucial for the 
ground state properties of the gas. However, when the gas is 
released from the trap, the radial expansion will be governed 
not by the interaction, as happened for the many channel case 
(ji ^ hw±) considered in Sec. [Ill but by the zero-point en- 
ergy associated with radial motion. Thus, for a strictly one- 
dimensional system it is not possible to generate a large phase 
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imprint from small initial density modulation during the radial 
expansion. Therefore, in this section as a means of impress- 
ing a phase we employ the second possibility, namely a short 
potential pulse, as mentioned in Sec. H] Phase imprinting of 
this kind has been proposed 1 13J and used extensively to gen- 
erate vortices in two and three dimensional condensates. The 
mechanism works as follows: starting at time — r we apply 
to the system, in its ground state, a short potential pulse, of 
duration r and of a prescribed spatial profile 



V(z,t) 



-T <t<0 



(24) 



where can be a deterministic or a random function of z. If 
the time interval r is shorter than the characteristic times of 
the system, then at the time instant t ~ the axial part of the 
many-body wavefunction will acquire a phase 



^'(2:1,. . . 2Ar;t=0) = exp 



N 



^oizi, ■ • • zn) 



(25) 

where <i>o(zi, . . . z^v) is the ground state wavefunction, prior 
to the action of the pulse, and is normalised to unity. (The 
complete wavefunction for the system in three dimensions of 
course also includes the radial factor Y[f=i XoiPj)-) At time 
t = 0, just after this phase has been impressed, the trapping 
potential is switched off and the gas undergoes radial expan- 
sion. The initial Gaussian function, Xo{p)^ will spread with 
time, retaining its Gaussian shape, and the density of the gas 
will evolve accordingly. Therefore we shall neglect interac- 
tions during the expansion, so that the gas is assumed to un- 
dergo free evolution and the z-dependent part of the many- 
body wavefunction evolves according to 



.^ d-^{zi,...ZN;t) 

at 



2ra dz: 



*(zi,...z7v;t) 



(26) 



which is to be solved with the initial condition given in 
Eq. dZSl l. This initial function contains all the information on 
the interacting groundstate, prior to the expansion. 

We are interested in the one-dimensional, z-dependent part 
of the particle density 



ni{z,t) = N 



.ZM]t)\ dz2 . 



N , , 
^ = J^F(z,t) 



(27) 

The actual, three-dimensional density n(p, z, t) is obtained 
by multiplying ni{z,t) by the radial factor |x(/c, t)P - the 
spreading Gaussian. Initially, F{z,t=0) — 1 but, as the ex- 
pansion proceeds, F{z,t) develops modulations in z, due to 
the initially impressed phase. In second quantised form 



L 



Fiz,t)^-m\z,tMz,t)\^) 



N 



(28) 



where I^I^) is the initial state vector, defined in position repre- 
sentation in Eq. dZSl l, while and ipizjt) are the free 



field operators 



-t-ikz i 



(29) 



with creation and annihilation operators aj, and that satisfy 

the commutation relation [a^, aj,] = dk,q. 
Using Eq. ( |29l ) and the relation 



(30) 

Eq. (|28] | can be written as 



F{z,t)^ J ^\Gp{z,t)\Mp) , (31) 
where n{p) is the momentum distribution function 

n{p) = ^ I dz{4>Hz)^ (0))e-'P^ (32) 



normalised so that 



(33) 



The function Gp{z, t) is 



Gp{z, t) — 



2Triht 



dC exp 



2hV ^' 



(34) 

In the mean field approach n{p) = 2'n5{p) and the expression 
for Gp{z,t) reduces to that for ip{z,t) given in Eq. (fTOl i of 
Sec. ini Such behaviour corresponds to the non-interacting 
limit of the Leib-Liniger gas which, under radial expansion, 
will exhibit large density variations and caustics, as discussed 
in Sec. Below we show that interactions inhibit caustic 
formation, and derive a condition for caustics to survive in the 
presence of interactions. 

The treatment of the integral in Eq. ( l34b is along the same 
lines as in Sec [III Caustics originate from rays emerging 
from points at which the second derivative of the phase 
in Eq. ( |34] | vanishes. Since this derivative does not con- 
tain p, we have the familiar condition for the singular points, 
d'^9{C) = —m/ht at C — C„. The first derivative, however, 
now contains p so that the analogue of Eq. (T3[ . for an arbi- 
trary function 6{(^) and keeping all physical dimensions, is 



Z = C+-(W)-P) 
m 



(35) 



This is essentially the definition of a ray, emerging from a 
point (, and having a particular value of p. The dependence 
on p implies that the rays emerging from the same singular 
point (n, but having different values of p, will not arrive at the 
point z at the same time. Thus focussing, which is the essence 
of the phenomenon of caustics, will be suppressed and caus- 
tics will get washed out. The condition for the existence of 
caustics follows from the integral in Eq. ( [34l i. Assuming that 
the function 9{C,) has a characteristic amplitude Oq and scale 
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of variation Rq, and expanding the phase in Eq. ( |34] | near a 
singular point (n, one obtains a contribution CO^R^'^rf' — prj 
where rj = C, — Cn C is a constant of order 1. Caustics 
originate from the rj'^-term and the relevant range of integra- 
tion is 77 < Ro9q . Therefore, if the contributing values 
of p, defined by the range over which n{p) is significant, are 

— 1/3 

such that RapOf^ <C 1, then caustics will survive. Since 
most of the weight of the momentum distribution n{p) is con- 
centrated in the interval p < 1/^ ifTil [Tsll . where ^ is the 
coherence length for the Bose gas, we arrive at the condition 

1 /3 

^0 ^ 1- It should be supplemented by the require- 

ment that, for a caustic to be formed and visible, there must 
be many particles in the region from which it originates. This 

— 1/3 

condition is Til i?o^'o ^ 1- Both conditions can be easily 
satisfied in the weak interaction limit where ^ni ^ 1. The sit- 
uation is different in the opposite case of strong interactions, 
or hard-core bosons. In this case ^ w l/j^i, and so caustics 
are absent in the strongly interacting limit. 



IV. TWO DIMENSIONAL SYSTEMS 

In this section we consider expansion of a quasi-two di- 
mensional condensate after release from a trap with strong 
axial and weak radial confinement, characterised as in SecHIl 
by frequencies uj^ and cjj^ but here with lOz ^ lo^. Again 
we neglect the weaker confinement, restricting our discus- 
sion to times t <^ l/wj, after release of the trap, and treat 
a BEC with initial density modulations, which are converted 
during expansion into an impressed phase. A difference be- 
tween the two-dimensional and one-dimensional geometries 
is that caustics in two dimensions consist of lines rather than 
points. Another difference is that, within the Thomas-Fermi 
approximation, a two-dimensional condensate in a disordered 
potential undergoes a percolation transition at a finite critical 
value of disorder strength. Because caustic formation is es- 
sentially a local pheomenon, the network of caustic lines does 
not show any critical behaviour that reflects this percolation 
transition, although it does change with disorder strength. 

The distinction between the initial and late phases of expan- 
sion is not as sharp in two-dimensions as it is in one dimen- 
sion. The reason for this is that the characteristic density near 
the center of the trap decreases with time t as in two di- 
mensions and as t^"^ in one dimension, with the result that the 
impressed phase grows logarithmically at long times in two 
dimensions, but reaches a limiting value in one dimension. 
We neglect this logarithmic growth and take the impressed 
phase to have a definite value 



(36) 



at times large compared to l/cj^, with characteristic amplitude 
6*0 and length scale 1 /fco. In this approximation the wavefunc- 
tion during the second phase of expansion can be factorised as 



For a two-dimensional system the axial part is given by a scal- 
ing function and our interest is in the evolution of the pla- 
nar part, $(rj^, t). In analogy with Eq. ( fTol i. it is given at late 
times by 



exp 



im 



i0(ri; 



(38) 

As in quasi-one dimensional systems, caustics are formed 
for 6*0 ^ 1 at values of the scaled time i > 1, and in this 
regime Eq. dSST l can be evaluated using the stationary phase 
method. The saddle points in r'^ are the solutions to 



ht 







(39) 



In the leading approximation, one such saddle point, at 
r'^ = r^, makes a contribution to $(rj^,t) of modulus 
[detAf(r'[)]-i/^ where 



M(rl) 



ht 



dy'dx'Oir'^, 



(40) 



Caustics stem from those saddle points at which detAf (r^) = 
0. Atoms in the expanding condensate originating from these 
points are focussed in such a way that the density |$(rj^, t)p 
is divergent within this approximation, which is equivalent to 
geometrical optics. To find the density in the vicinity of caus- 
tics, it would be necessary to take into account higher deriva- 
tives of 0{t±^ when calculating the integral in Eq. ( [38] ). We 
do not do this, restricting ourselves instead to a discussion of 



the positions of caustics. The condition detAf(rj 



de- 



(37) 



fines a set of lines in the atomic cloud after the initial phase of 
expansion, which give rise to caustics. Both the number and 
the shape of these lines in the initial plane depend on the time 
at which density in the expanding cloud is to be measured, 
but they have fixed limits for t 3> 1. Atoms starting from 
points on these lines travel with velocity {h/m)V9, and reach 
points in the cloud at time t with coordinates given by the 
solutions to Eq. ( |39] l. In this way lines in the initial plane 
are mapped to moving lines of high density in the expanding 
cloud. 

To illustrate these general ideas, consider the example of a 
periodic impressed phase 

9{v\^ = 6 q{cos kox + cos kf)y) . (41) 

In this case the condition detA/(r^) — yields 

{i-^ - cosfcoa;*)(f-i - cos hoy*) = , (42) 

defining two sets of parallel lines in the initial plane, x* — 
fcfj"^ cos~^(l/f) and y* = k^^ cos^^{l/i). Caustics derive 
from these, forming two similar sets of lines in the final plane, 
X = X* — t sin kox* and y ^ y* ~ t sin k^y* respectively. 

The consequences in two dimensions of a potential V{y±^ 
strong enough to generate significant density modulations be- 
fore expansion can be discussed for large 6*0 using the same 
approach as in Sec. [Ill In this way we find that two conditions 
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FIG. 5: Effect of decreasing /i in the condensate before expan- 
sion, on caustic formation in two dimensions, illustrated for a pe- 
riodic potential V{r±). Each of the four panels shows the ini- 
tial plane. Caustics evolve from lines with detM(rx) — and 
|$(rx,0)p > 0. These are marked in bold, while lines with 
detA/(rx) = and |<&(rx,0)|^ = are shown dotted. Panel (a): 
large jj., |<I>(rx, 0)]'^ > everywhere; (b) on reducing /i, the con- 
densate develops small holes but percolates, while caustics become 
fragmented; (c) reducing /i further, the condensate reaches the perco- 
lation threshold, but caustics show no signature of this; (d) reducing 
fi still further, caustics shrink and eventually disappear. 



function of position and time, the effects of particle statistics 
enter the main result, Eq. dSTT i. only through the momentum 
distribution of particles before expansion. The criterion for 
formation of caustics in the expanding Fermi gas is therefore 
the same as for the Lieb-Liniger gas, but with the Fermi wave- 
length Ap taking the place of the coherence length ^ in the 
expressions given in Sec. |III] Hence caustics are absent from 
the Fermi gas, for the same reason as in the Bose gas with 
hard-core interactions. There are nevertheless differences be- 
tween non-interacting fermions and hard-core bosons. They 
stem from the Fermi surface discontinuity in the momentum 
distribution. Within the approximations of geometrical optics, 
this discontinuity leads to sharp peaks in the derivative of the 
relative density F(z, t) with respect to position z or time t. To 
show this, we note from Eq. ( [34] i that 



d,\Gpiz,t)\^ = -dp\Gpiz,t)\^ . (43) 
With Eq. (EB this yields 

d.Fiz,t) = -^J ^\G,{z,t)\'d,Mp) 

= ^[|Gp,(z,i)p-|G_,,(z,i)p] (44) 

where pp is the Fermi wavevector The behaviour of Gp{z, t) 
has been analysed in Sec. Ull the value of p influences the po- 
sition of caustics but not their formation. The Fermi gas there- 
fore shows the same extrema in the derivative of the density 
as are found for a BEC in the density itself. 

VI. SUMMARY 



must be satisfied in order that a line segment in the initial plane 
will give rise to a caustic after expansion: it is necessary, first, 
that detM(rj^) = and, second, that the initial relative den- 
sity |$(rj^,0)p is non-zero on the line. In consequence, as 
the potential strength is increased, or the chemical potential is 
reduced, caustic lines first develop breaks, and then disappear 
altogether Such an evolution with decreasing fi of the lines in 
the inital plane that generate caustics is illustrated in Fig|5]for 
the periodic potential that underlies Eq. dTIT l. 

The theory of caustic formation resulting from a ran- 
dom phase 6'(rx) is analogous to the treatment of the phase 
screen problem, which has been studied extensively for two- 
dimensional systems in the context of optics [2]. In particular, 
the morphology of caustic lines for the random case is dis- 
cussed in lfl6ll . 



V. FERMIONS 

It is interesting to ask about problems similar to the ones 
we have discussed, but with fermions in place of bosons. To 
be specific, consider expansion of a strictly one dimensional 
system with an imprinted phase, as in Sec. |III1 but for non- 
interacting fermions rather than interacting bosons. Restrict- 
ing our attention to the expectation value of the density as a 



We have discussed free motion of one and two dimensional 
atomic gases with an initial impressed phase that varies peri- 
odically or randomly as a function of position. Gradients of 
this phase represent initial velocities and lead to density varia- 
tions that grow with time. The characteristic amplitude and 
size Rq of spatial variations of this phase are key parameters. 
The limit 0o ^ 1 is both the most interesting regime, because 
density maxima are largest, and a tractable one theoretically, 
because a treatment analogous to geometrical optics provides 
the leading approximation. The evolution of atomic density 
fluctuations with time has close links to problems in optics 
involving caustic formation. In the context of atomic gases, 
caustics are maxima of density, near points in one dimensional 
systems or along lines in the two dimensional case. For atoms 
of mass m they form at a characteristic time t* — mR^JhOf) 
and at longer times t the density within a caustic decays as 
t^^. Since caustics originate from small regions of the initial 
atomic cloud, variations in the initial density simply modu- 
late the density on caustics. In particular, caustic formation 
is suppressed during expansion of a fragmented condensate if 
the initial density is zero at points that would otherwise be the 
origin for caustics. 

We have argued that a recent experiment IH] in which large 
density modulations are observed in an elongated BEC after 
release from a disordered potential should be understood in 
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terms of caustic formation. For the future it would be of in- SRC under Grant No. EP/D050752/1. 
terest to design experiments with larger values of 6*0 for both 
one and two dimensional systems. 
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